{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.1.3 Reissner-Mindlin plate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To avoid shear locking when the thickness $t$ becomes small several methods and elements have been proposed. We discuss two approaches:\n",
    "\n",
    "1: Mixed Interpolation of Tensorial Components (MITC)\n",
    "\n",
    "2: Rotations in Nèdèlec space using the Tangential-Displacement Normal-Normal-Stress method (TDNNS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixed Interpolation of Tensorial Components (MITC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.geom2d import SplineGeometry\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up benchmark example with exact solution. Young modulus $E$, Poisson ratio $\\nu$, shear correction factor $\\kappa$ and corresponding Lamè parameters $\\mu$ and $\\lambda$. Geometry parameters are given by thickness $t$ and radius $R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu, k = 10.92, 0.3, 5/6\n",
    "mu  = E/(2*(1+nu))\n",
    "lam = E*nu/(1-nu**2)\n",
    "#shearing modulus\n",
    "G = E/(2*(1+nu))\n",
    "#thickness, shear locking with t=0.1\n",
    "t = 0.1\n",
    "R = 5\n",
    "#force\n",
    "fz = 1\n",
    "\n",
    "order = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to symmetry we only need to mesh one quarter of the circle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sg = SplineGeometry()\n",
    "pnts = [ (0,0), (R,0), (R,R), (0,R) ]\n",
    "pind = [ sg.AppendPoint(*pnt) for pnt in pnts ]\n",
    "sg.Append(['line',pind[0],pind[1]], leftdomain=1, rightdomain=0, bc=\"bottom\")\n",
    "sg.Append(['spline3',pind[1],pind[2],pind[3]], leftdomain=1, rightdomain=0, bc=\"circ\")\n",
    "sg.Append(['line',pind[3],pind[0]], leftdomain=1, rightdomain=0, bc=\"left\")\n",
    "mesh = Mesh(sg.GenerateMesh(maxh=R/3))\n",
    "mesh.Curve(order)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on the boundary conditions (simply supported or clamped) we have different exact solutions for the vertical displacement. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = sqrt(x**2+y**2)\n",
    "xi = r/R\n",
    "Db = E*t**3/(12*(1-nu**2))\n",
    "\n",
    "clamped = True\n",
    "\n",
    "#exact solution for simply supported bc\n",
    "w_s_ex = -fz*R**4/(64*Db)*(1-xi**2)*( (6+2*nu)/(1+nu) - (1+xi**2) + 8*(t/R)**2/(3*k*(1-nu)))\n",
    "#Draw(w_s_ex, mesh, \"w_s_ex\")\n",
    "#exact solution for clamped bc\n",
    "w_c_ex = -fz*R**4/(64*Db)*(1-xi**2)*( (1-xi**2) + 8*(t/R)**2/(3*k*(1-nu)))\n",
    "Draw(w_c_ex, mesh, \"w_c_ex\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use (lowest order) Lagrangian finite elements for rotations and the vertical deflection together with additional internal bubbles for order >1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if clamped:\n",
    "    fesB = VectorH1(mesh, order=order, orderinner=order+1, dirichletx=\"circ|left\", dirichlety=\"circ|bottom\", autoupdate=True)\n",
    "else:\n",
    "    fesB = VectorH1(mesh, order=order, orderinner=order+1, dirichletx=\"left\", dirichlety=\"bottom\", autoupdate=True)\n",
    "fesW = H1(mesh, order=order, orderinner=order+1, dirichlet=\"circ\", autoupdate=True)\n",
    "fes = FESpace( [fesB, fesW], autoupdate=True ) \n",
    "(beta,u), (dbeta,du) = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Direct approach where shear locking may occur\n",
    "$$\\frac{t^3}{12}\\int_{\\Omega} 2\\mu\\, \\varepsilon(\\beta):\\varepsilon(\\delta\\beta) + \\lambda\\, \\text{div}(\\beta)\\text{div}(\\delta\\beta)\\,dx + t\\kappa\\,G\\int_{\\Omega}(\\nabla u-\\beta)\\cdot(\\nabla\\delta u-\\delta\\beta)\\,dx = \\int_{\\Omega} f\\,\\delta u\\,dx,\\qquad \\forall (\\delta u,\\delta\\beta). $$\n",
    "\n",
    "Adding interpolation (reduction) operator $\\boldsymbol{R}:[H^1_0(\\Omega)]^2\\to H(\\text{curl})$. Spaces are chosen according to [<a href=\"http://math.aalto.fi/~rstenber/Publications/M3AS91.pdf\">Brezzi, Fortin and Stenberg. Error analysis of mixed-interpolated elements for Reissner-Mindlin plates. <i>Mathematical Models and Methods in Applied Sciences 1</i>, 2\n",
    "  (1991), 125-151.</a>]\n",
    "$$\\frac{t^3}{12}\\int_{\\Omega} 2\\mu\\, \\varepsilon(\\beta):\\varepsilon(\\delta\\beta) + \\lambda\\, \\text{div}(\\beta)\\text{div}(\\delta\\beta)\\,dx + t\\kappa\\,G\\int_{\\Omega}\\boldsymbol{R}(\\nabla u-\\beta)\\cdot\\boldsymbol{R}(\\nabla\\delta u-\\delta\\beta)\\,dx = \\int_{\\Omega} f\\,\\delta u\\,dx,\\qquad \\forall (\\delta u,\\delta\\beta). $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gamma = HCurl(mesh,order=order,orderedge=order-1, autoupdate=True)\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += t**3/12*(2*mu*InnerProduct(Sym(grad(beta)),Sym(grad(dbeta))) + lam*div(beta)*div(dbeta))*dx\n",
    "#a += t*k*G*InnerProduct( grad(u)-beta, grad(du)-dbeta )*dx\n",
    "a += t*k*G*InnerProduct(Interpolate(grad(u)-beta,Gamma), Interpolate(grad(du)-dbeta,Gamma))*dx\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += -fz*du*dx\n",
    "\n",
    "gfsol = GridFunction(fes, autoupdate=True)\n",
    "gfbeta, gfu = gfsol.components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define function for solving the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP():\n",
    "    fes.Update()\n",
    "    gfsol.Update()\n",
    "    with TaskManager():\n",
    "        a.Assemble()\n",
    "        f.Assemble()\n",
    "        inv = a.mat.Inverse(fes.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "        gfsol.vec.data = inv * f.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Solve, compute error, refine, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "l = []\n",
    "for i in range(5):\n",
    "    print(\"i = \", i)\n",
    "    SolveBVP()\n",
    "    if clamped:\n",
    "        norm_w = sqrt(Integrate(w_c_ex*w_c_ex, mesh))\n",
    "        err = sqrt(Integrate((gfu-w_c_ex)*(gfu-w_c_ex), mesh))/norm_w\n",
    "    else:\n",
    "        norm_w = sqrt(Integrate(w_s_ex*w_s_ex, mesh))\n",
    "        err = sqrt(Integrate((gfu-w_s_ex)*(gfu-w_s_ex), mesh))/norm_w\n",
    "    print(\"err = \", err)\n",
    "    l.append ( (fes.ndof, err ))\n",
    "    mesh.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convergence plot with matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"relative error\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "plt.ion()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TDNNS method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Invert material law of Hooke and reset mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CMatInv(mat, E, nu):\n",
    "    return (1+nu)/E*(mat-nu/(nu+1)*Trace(mat)*Id(2))\n",
    "\n",
    "mesh = Mesh(sg.GenerateMesh(maxh=R/3))\n",
    "mesh.Curve(order)\n",
    "Draw(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of using Lagrangian elements together with an interpolation operator we can directly use H(curl) for the roation $\\beta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=1\n",
    "if clamped:\n",
    "    fesB = HCurl(mesh, order=order-1, dirichlet=\"circ\", autoupdate=True)\n",
    "    fesS = HDivDiv(mesh, order=order-1, dirichlet=\"\", autoupdate=True)  \n",
    "else:\n",
    "    fesB = HCurl(mesh, order=order-1)\n",
    "    fesS = HDivDiv(mesh, order=order-1, dirichlet=\"circ\", autoupdate=True)\n",
    "fesW = H1(mesh, order=order, dirichlet=\"circ\", autoupdate=True)\n",
    "\n",
    "fes = FESpace( [fesW, fesB, fesS], autoupdate=True ) \n",
    "(u,beta,sigma), (du,dbeta,dsigma) = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the TDNNS method with the stress (bending) tensor $\\boldsymbol{\\sigma}$ as proposed in [<a href=\"https://link.springer.com/article/10.1007/s00211-017-0883-9\">Pechstein and Schöberl. The TDNNS method for Reissner-Mindlin plates. <i> Numerische Mathematik 137</i>, 3 (2017), 713-740</a>]\n",
    "\\begin{align*}\n",
    "&\\mathcal{L}(u,\\beta,\\boldsymbol{\\sigma})=-\\frac{6}{t^3}\\|\\boldsymbol{\\sigma}\\|^2_{\\mathcal{C}^{-1}} + \\langle \\boldsymbol{\\sigma}, \\nabla \\beta\\rangle + \\frac{tkG}{2}\\|\\nabla u-\\beta\\|^2 -\\int_{\\Omega} f\\cdot u\\,dx\\\\\n",
    "&\\langle \\boldsymbol{\\sigma}, \\nabla \\beta\\rangle:= \\sum_{T\\in\\mathcal{T}}\\int_T\\boldsymbol{\\sigma}:\\nabla \\beta\\,dx -\\int_{\\partial T}\\boldsymbol{\\sigma}_{nn}\\beta_n\\,ds\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = specialcf.normal(2)\n",
    "    \n",
    "a = BilinearForm(fes)\n",
    "a += (-12/t**3*InnerProduct(CMatInv(sigma, E, nu),dsigma) + InnerProduct(dsigma,grad(beta)) + InnerProduct(sigma,grad(dbeta)))*dx\n",
    "a += ( -((sigma*n)*n)*(dbeta*n) - ((dsigma*n)*n)*(beta*n) )*dx(element_boundary=True)\n",
    "a += t*k*G*InnerProduct( grad(u)-beta, grad(du)-dbeta )*dx\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += -fz*du*dx\n",
    "\n",
    "gfsol = GridFunction(fes, autoupdate=True)\n",
    "gfu, gfbeta, gfsigma = gfsol.components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP():\n",
    "    fes.Update()\n",
    "    gfsol.Update()\n",
    "    with TaskManager():\n",
    "        a.Assemble()\n",
    "        f.Assemble()\n",
    "        inv = a.mat.Inverse(fes.FreeDofs())\n",
    "        gfsol.vec.data = inv * f.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve, compute error, and refine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []\n",
    "for i in range(5):\n",
    "    print(\"i = \", i)\n",
    "    SolveBVP()\n",
    "    if clamped:\n",
    "        norm_w = sqrt(Integrate(w_c_ex*w_c_ex, mesh))\n",
    "        err = sqrt(Integrate((gfu-w_c_ex)*(gfu-w_c_ex), mesh))/norm_w\n",
    "    else:\n",
    "        norm_w = sqrt(Integrate(w_s_ex*w_s_ex, mesh))\n",
    "        err = sqrt(Integrate((gfu-w_s_ex)*(gfu-w_s_ex), mesh))/norm_w\n",
    "    print(\"err = \", err)\n",
    "    l.append ( (fes.ndof, err ))\n",
    "    mesh.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convergence plot with matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"relative error\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.ion()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
